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Abstract 

Two-dimensional and three-dimensional kinetic simulation results reveal the importance of 
the Lower-Hybrid Drift Instability LHDI to the onset of magnetic reconnection. Both explicit 
and implicit kinetic simulations show that the LHDI heats electrons anisotropically and increases 
the peak current density. Linear theory predicts these modifications can increase the growth 
rate of the tearing instability by almost two orders of magnitude and shift the fastest growing 
modes to significantly shorter wavelengths. These predictions are confirmed by nonlinear kinetic 
simulations in which the growth and coalescence of small scale magnetic islands leads to a rapid 
onset of large scale reconnection. 
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I. INTRODUCTION 



Understanding the onset of magnetic reconnection remains an open unsolved question. 
Reconnection is observed in a great variety of systems, from special purpose laboratory 
experiments 

mm 

, to fusion devices [3] and space p, lD, |8| and astrophysical plasmas 
BHQ. Ou. focus is on .econnection in the n,agneto.phe.e due to u.te.actiou. with the 
solar wind. Reconnection develops mainly in current sheets at the magnetopause j7( and in 
the tail Q. 

In both situations, the field can be approximated by a Harris current sheet In fact, a 
reference configuration for such plasmas is proposed, the so-called GEM challenge, to b ring 
uniformity, repeatability and ease of inter-comparison among various models (see Ref. 13| 
and references therein). 

The original GEM challenge is initialized with a large magnetic island (perturbation of 
10% amplitude) so that both magnetohydrodynamic (MHD) modeling and kinetic simula- 
tions can study reconnection, after reconnection starts. Thus, the GEM challenge bypasses 
reconnection onset and allows researchers to study the physics of reconnection in its nonlinear 
phase. In our simulations, we consider a GEM equilibrium in which an initial perturbation 
is absent so that we can study reconnection onset. 

There have been many previous studies of reconnection onset. Among kinetic studies, 
which are most relevant to coUisionless magnetospheric plasmas where the current sheet 
thicknesses are of the order of the ion gyroradius or ion skin depth, the onset of reconnection 
has been long attributed to the tearing instability j^, 15, 16 1. 



However, the tearing instability saturates at low levels: For typical magnetospheric con- 
figurations, tearing instability saturation levels are so low that they yield tiny magnetic 
islands with widths of the order of the electron skin depth |ll^. These are far too small 
to trigger reconnection by any of the recently investigated, nonlinear mechanisms for fast 
reconnection The GEM challenge, with its large initial perturbation, produces a large 
enough starting island size to decouple electron and ion motion. This allows fast reconnec- 
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23|. However, without 



tion through Hall or whistler physics reconnection 
an initial perturbation, it is difficult for tearing modes alone to produce sufficiently large 
islands, and fast reconnection simply does not occur. Conceivably, the tiny islands could 
coalesce into larger islands until they reach the critical size necessary for fast reconnection 



physics. However, this mechanism is quite slow, and the present paper seeks other mecha- 
nisms that occur on much faster time scales. The tearing instability is also believed to be 
suppressed by the presence of even a relatively small vertical magnetic field (i.e., orthogonal 
to the current sheet) 
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28| . We defer this issue to a later study. 
Our study of the onset of coUisionless magnetic reconnection in a GEM challenge Harris 
sheet equilibrium without an initial perturbation shows that the Lower Hybrid Drift 
Instability (LHDI) (see Ref. 29| for a recent review of the literature) has an important role 
in the onset of reconnection in three dimensions. In two dimensions, tearing instabilities 
saturate at low amplitudes, as expected. In three dimensions, the LHDI modifies the current 
sheet in ways that cause the tearing mode to grow more vigorously, as also observed by Sholer 



et al. [30[, and by Shinohara and Fujimoto |31|. 

Waves at the lower hybrid frequency are observed near reconnection sites 
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and 



are often considered as a source of anomalous resistivity. In the earliest model of reconnection 
presented by Sweet and Parker, based on resistive MHD, the classical Spitzer resistivity is 
not able to explain the reconnection rates observed. Microinstabilities can enhance the 
plasma resistivity beyond the classical values, and thus, the LHDI, which is unstable to 
a broad range of wavelengths and frequencies, has been studied extensively as a source of 
anomalous resistivity. The fastest growing modes have wavelengths kyPe ~ 1, ~ il/, 
are active in the low-/5 region of the current sheet and are predominantly electrostatic 
Longer wavelength modes {ky^fpipl ^ 1, 7 ^ Vtci) have a large electromagnetic cornponent 



and for sufficiently thin sheets can penetrate into the center of the current sheet [36|, Is? ]. 
However, observations show that the LHDI is confined to the edge of current sheets, its 
amplitude seems uncorrelated with the reconnection electric field, and it saturates at too 
ow amplitude to explain the enhanced reconnection rate through an anomalous resistivity 
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Recently, the effects of the Lower-Hybrid drift instability have been revisited. It has 
been observed that in the non-linear phase of the evolution of the LHDI in a Harris sheet, 
modifications of the profiles are induced. 



First, the current sheet is thinned and the electron current profile is peaked 
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The effect of electron current shee 
to electron density modification 



rgeaking is primarily due to electron acceleration and not 
From these results the possibility emerges that the 
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modifications of tlie current profile induced by the LHDI allow the onset of secondary insta- 
bilities and modify the growth rates of instabilities already present in the initial equilibrium. 

Second, previous results based on two-dimensional simulations in the plane of the LHDI 
(orthogonal to the magnetic field and to the plane where tearing develops) have indeed shown 
that a Kelvin-Helmholtz Instability (KHI) develops at much faster rates than predicted by 



linear theories based on the initial equilibrium (see Ref. 



□ 



and references therein). 
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Furthermore, previous three-dimensional simulations 
modification of the initial profile caused by the LHDI can allow the onset of reconnection 
and enhance the growth of the tearing instability. 

In the present paper we revisit the issue with more detailed simulations and with more 
extensive diagnostics, by which we uncover new physics and confirm previous results. 

First, we discover that the nonlinear evolution of the LHDI not only heats electrons (an 
effect documented in many previous works, e.g. Ref. js^) but preferentially heats electrons 



perpendicular to the magnetic fie 
of reconnection. Previous work 



d. Anisotropic heating has a great impact on the onset 



45 



46 



suggests that systems with an electron 



(or an ion) anisotropy are more unstable to the tearing instability when the perpendicular 
temperature exceeds the parallel temperature. We conduct a number of simulations and 
we use the linear Vlasov theory to calculate the effect of temperature anisotropy on tearing 
onset. 

Second, we confirm that the current sheet is thinned by the LHDI. The thinning can 
be explained by nonlinear effects of the LHDI [48 1. The thinning has the direct effect of 
promoting the growth of the tearing instability. We conduct a number of simulations and 
theoretical studies based on linear theory to estimate the effects of current sheet thinning 
on the linear growth rate. 

In Sec. II, we describe our physical system and numerical approach. In Sec. HI, we 
present our simulation results. In Sec. IV, we discuss our results and conclusions. 
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II. PHYSICAL SYSTEM AND THE NUMERICAL APPROACH 



A Harris current sheet equilibrium is considered 1^ , with an initial magnetic field given 

by 

Bo(2;) = 5otanh(z/L)e, (1) 

and a plasma density given by 

no{z) = riQ secli^ (z / L) (2) 

In the present paper, similar parameters to the GEM challenge are used [l^. In the 
GEM challenge, the current sheet thickness is L = O.Srfj, the temperature ratio is Tj/Tg = 5, 
and the ion drift velocity in the y direction is V^o = 1-67Va. Rather than the standard mass 
ratio, mi/rrie = 180 is used. The Alfven velocity, Va, and the ion inertial length, di = c/upi, 
are defined with the density no and the field Bq. In contrast to the GEM challenge, no 
background density is introduced. This is an important point, since the presence of a 
background is strongly stabilizing to the LHDI. 

The standard dimensionless parameters necessary to characterize a Harris current sheet 
can thus be summarized as, Upe/Qce = 2.88, pi/L = 1.828, mi/me = 180, and Ti/T^ = 5, 
with vth,s = \/2Tjms. 

The boundary conditions for the particles and fields are periodic in the x and y directions. 
Conducting boundary conditions are imposed for the fields at the z boundaries while refiect- 
ing boundary conditions are used for the particles. In order to study the reconnection onset 
and in contrast to the GEM challenge jl3], the Harris equilibrium is initially unperturbed. 

We simulate the dynamics in both the reconnection, {x,z), and current-aligned, {y,z), 
planes. In the {x, z) plane, the domain is x = 25.6L x 12. 8L, corresponding to 
Lx X Lz = I2.8di X 6.4dj. (The GEM challenge domain is L^x Lz = 51.2L x 25. 6L.) In the 
{y, z) plane, the domain is Ly x = 32L x 12.8L. We also perform a three-dimensional 
simulation, for which the domain is x Ly x = 25.6L x 16L x 12. 8L. (This simulation 
is performed with CELESTE3D only, which is described below.) 

To investigate the evolution of the system, a kinetic linear code and two nonlinear PIC 
simulation codes are used. The linear Vlasov results for the tearing tearing mode are cal- 
culated us,,. t>,e fo™,..^ exact app.oacK desc.bed ,„ Q. TKis tech.,ue employs a 
normal mode calculation using a full Vlasov description for both ions and electrons. The 



orbit integrals arising from the linear Vlasov theory are treated numerically using the exact 
unperturbed particle orbits and including the form of the perturbation inside the integral. 
Both electromagnetic and electrostatic contributions to the field perturbation are retained 
and resulting system of integro-differential equations is solved using a finite element expan- 
sion of the eigenfunction j37|. The basic strategy involves a normal mode calculation for 
perturbations of the form 

= 0(z) exp{—iut + ikyH + ik^x) , (3) 
A = A{z) exp{—iujt + ikyH + ik^x) , 

where the complex functions 0, A are the perturbed electrostatic and vector potentials. For 
a given Vlasov equilibrium and for a given choice of wavector {kx,ky), the code computes the 
real frequency, growth rate (real and imaginary part of u) and the complex eigenfunctions 
(j){z) and A{z) which describe the mode structure. 

The nonlinear dynamics are simulated by two PIC codes, an explicit simulation code 
NPIC, and an implicit simulation code CELESTE3D. 

The explicit plasma simulation code NPIC is based on a well-known explicit electro- 



magnetic algorithm 



5G 



5l|. The particle trajectories within NPIC are advanced using the 
leapfrog technique, and particle moments are accumulated with area weighting. The simu- 
lations are run on a parallel computer, using domain decomposition with calls to the MPI 
library. (In the present work, the explicit simulations are run on the Los Alamos Q-machine 
using as many as 128 nodes.) 



The implicit plasma simulation code, CELESTE3D, solves the full set of Maxwell- Vlasov 
equations using the implicit moment method 0, 0, |54|. Both Maxwell's and Newton's 
equations are discretized implicitly in time. The implicit simulations are run on a worksta- 
tion. 

The non-linear simulations are performed by the two codes with very different simulation 
parameters. In the tearing plane, NPIC employs a x = 1280 x 640 grid, a time step 
OceAt = 0.03, and 160 ■ lO'' particles. CELESTE3D uses a A^^ x A^^ = 64 x 64 grid, with 
time step ^Ice^t = 0.45, and a total of 5 ■ 10^ computational particles. In the current aligned 
plane, NPIC employs Ny x = 1600 x 640 grid, a time step Qce^t = 0.03, and 2 ■ 10® 
particles , while CELESTE3D uses a A^^. x = 128 x 64 grid, with time step Qce^t = 0.7, 
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and a total of 1 ■ 10^ computational particles. The simulation parameters are summarized 
in Table I. 

Detailed comparison of the plasma dynamics in NPIC and CELESTE3D are made in a 
study of the non-linear phase of magnetic reconnection in plasmas with different (3 values 
2^. The comparison shows that the physical mechanisms revealed by the two codes agree, 
which increases confidence in their validity. Similar comparisons are made with other codes 
for the GEM challenge, which show that results with CELESTE3D are comparable in detail 

0|. The same kinds of comparisons are made in this 



with those of explicit simulations |54 
study. 



III. SIMULATION RESULTS 

To study the onset of reconnection, we simulate the reference configuration defined above, 
similar to the GEM challenge 13], but without initial perturbation, without background 
plasma, mj/mg = 180, and smaller box size. It is left to the natural noise of PIC to excite 
any instability. The result is strikingly different from published GEM challenge results. 
The tearing mode saturates at low amplitude, and the reconnected flux is a small fraction 
of the available. If a third dimension in the current aligned direction, is added to the 
simulation of the same physical system, the tearing mode grows to much larger amplitude 
before saturating and reconnection occurs. Two effects of the LHDI instability appear to 
cause a dramatic change in both the linear and the non-linear phase of the tearing instability 
in three-dimensional dynamics, compared with the two-dimensional dynamics. They are 
anisotropic electron heating and current sheet thinning. Both of these effects enhance the 
linear growth rate of the tearing instability and strongly affect and increase the saturation 
amplitude. In the following subsections, the results of the simulations are described in detail. 

A. Linear tearing and saturation 

1. Two-dimensional simulation 

We consider a two-dimensional simulation in the tearing plane (x, z), with the parameters 
defined above, similar to the GEM challenge. For this configuration, linear theory predicts 
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that the fastest growing tearing mode has wavelength kj.L = 0.5 and growth rate 7 = 
0.176Qci- In Fig- 1, the amphtude of tearing modes with rrix = 1, 2, 3, and 4, corresponding 
to kr^L = 0.25, 0.5, 0.75, and 1, are shown as a function of time. The growth of higher mode 
numbers is neghgible. 

At the beginning of the simulation. Fig. 1 shows the fast growth of the tearing mode 
with = 2, k^L = 0.5, in agreement with the linear theory. The growth of the mode with 
TUx = 1 reveals the merging of two islands. 

The tearing instability saturates at a low level. The half-width of the island at saturation 
is w ~ 0.46L at time tQ^i = 83. By comparison, GEM challenge reconnection encompasses 
the whole domain at tQ^i ~ 30 (w ~ lOL). CELESTE3D and NPIC agree in showing low 
amplitude saturation of the tearing mode, but the coarser grid spacing used in CELESTE3D 
does not allow an accurate estimate of the island half width. It should be noted that the 
growth of tearing modes in NPIC with rrix = 3, 4 at tQci ~ 80 is possibly due to the numerical 
heating always present in explicit PIC simulations like NPIC, which can in principle affect 
the results of long time scale simulations. In any case, this slow growth at late time in the 
two-dimensional tearing simulation is several orders of magnitude slower than the growth 
and non-linear development of the LHDI in three-dimension that is the subject of the present 
work. 

The saturation of tearing has been studied theoretically in high and low-jS plasmas (e.g., 
see Refs. 0,y,0)- 

In high-/? plasmas, like the ones considered in the GEM challenge, 
anisotropic heating of electrons causes saturation of tearing ^^l- During the growth of 
the tearing instability, (Teyy + T^zz) I'^T^xx (defined as Tej_/Te||) is reduced below 1, and 
Te^±/Te\\ < 1 has a strongly stabilizing effect on tearing. In fact, in the two-dimensional 
simulations of tearing at time tVLd = 80, NPIC gives Te±/Te\\ = 0.87 and CELESTE3D gives 
Te±/Te\\ = 0.83. For Te±/Te\\ = 0.87 the linear code predicts the mode 171^ = 2 {k^L = 0.5) 
to be stable and the mode rrix = 1 {kxL = 0.25) to be weakly unstable (7 = O.OAQci)- 
However, this growth rate is based on a bi-Maxwellian electron distribution. The simulated 
velocity distribution is more complex and might yield a different growth rate. 

It should be remarked that in simulations performed in a more extended domain, that 

n 

small islands can coalesce repeatedly into larger islands [58]. If the system is big enough, 
the final merged island can grow enough to allow the fast reconnection physics shown by 
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the GEM challenge. However, this coalescence process occurs on a long time scale compared 
to the time scales typical of the three-dimensional instabilities that are the subject of the 
present manuscript. 



2. Three-dimensional simulation 

CELESTE3D is used to perform a three-dimensional simulation of the Harris current 
sheet described in Sect. H, in a domain with x Ly x = 25. 6L x 16L x 12. 8L. Tear- 
ing, current aligned (i.e., LHDI and KHI), and oblique instabilities are all involved in the 
evolution of the system. 

The results of three-dimensional simulations performed by CELESTE3D have been di- 
rectly compared with satellite observations j^. A close agreement has been observed be- 
tween simulation and satellite observations, concerning current sheet kinking, current bifur- 
cation, and reconnecting modes. 

In Fig. 2, the amplitude of the Fourier component of the magnetic field is plotted as a 
function of the mode number rrix (aligned with the magnetic field) and rriy (aligned with the 
current) at sequential times. The plots average over z and Fast Fourier Transform (FFT) 
in (x, y) . The average over z suppresses odd parity modes in and thus only even parity 
modes, such as KHI and LHDI, appear in Fig. 2. Since B^ and B^ have opposite parity, 
the even modes, like the tearing modes (not shown in Fig. 2) are picked up by z-averages 
of Fourier modes of Bz- 

Fig. 2 shows the development of the initial noise, which excites mostly low mode number 
instabilities m^^ and my at tQci = 0.2. The presence of the electrostatic LHDI is reflected at 
time tQci = 2.5 through a significant amplitude of the mode number = 0, rriy = 12 — 14, 
which corresponds to kyp^ ^ 0.2. The electrostatic LHDI, studies of which are summarized in 



Ref. |29|, is dominated by an electrostatic component on the edge of the current sheet. The 
characteristic wavelength is of order kyPe ^ 1, and frequency u ^ Qih The electrostatic 
LHDI instability causes a velocity shear that enhances the growth rate of the electromagnetic 
LHDI and causes a KHI to grow |4ll.l4^. Later, at tQci = 5, the maximum amplitude modes 
are m^, = and rUy = 5 — 6, which corresponds to ky^fpipl ^ 0.5). These longer wavelength 
modes are the electromagnetic LHDI, which grows at center of the current sheet and has 
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wavelength ky^JJ^pl ^ 1 and growth rate 7 ^ Vld [37]. Velocity shear created by the LHDI 
triggers a KHI at a still later time, so that by tQci = 20 the current sheet is dominated by 
a single, domain-sized kink with mode numbers = and rriy = 1 {kyL ^ 0.4). 

Figure 3 shows the amplitude of the Fourier modes of the field as a function of 
and rUy at different times. Again, averaging over z eliminates the odd parity modes and 
leaves only even parity modes. The modes shown in Fig. 3 are complementary to the those 
in Fig. 2. The Bz Fourier components are signatures of the tearing instability. The tearing 
instability starts to grow significantly at tQci = 5, with mode number = 4 and my = 
{kxL = 1, a higher frequency than the expected from the linear theory). Subsequently, the 
4 magnetic islands merge into 2 islands at tQci = 10, and then into a 1 that encompasses 
the whole domain by tQd ~ 20. This is comparable to the reconnection time shown by 
the GEM challenge simulation in two dimensions with a large initial perturbation jl^. No 
significant growth of oblique tearing modes is observed during the simulation. 

Two main conclusions can be drawn from the three-dimensional simulation. First, in 
contrast to simulations in two-dimensions, the tearing mode does not saturate at a low 
amplitude in three dimensions. Instead, it encompasses the whole domain in a time com- 
parable to the GEM challenge simulation with a large initial perturbation Second, no 
significant growth of oblique modes is seen in the simulation. 

Thus, it can be argued that the current sheet dynamics depends only on some nonlinear 
interaction between the tearing instability, which develops in the (x, z) plane, and current 
aligned instabilities, which grow in the {y, z) plane. As the time scale of instabilities in 
the {y, z) plane is an order of magnitude faster than the tearing instability, it follows that 
the fundamental dynamics of a three-dimensional current sheet can be understood more 
conveniently by performing simulations in the (x, z) and {y, z) plane, and analyzing the 
effects of the very rapid current aligned instability on the more slowly growing tearing 
mode. 

B. Current sheet modifications by the LHDI 

In simulations performed in the current aligned, {y,z), plane, two modifications of the 
initial equilibrium appear to enhance the tearing growth rate. Both are non-linear conse- 
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quences of the LHDI. 

The LHDI causes anisotropic heating of the electrons. In Fig. 4, the y-averaged electron 
temperature ratio, Te^/Te||, is plotted as a function of z at t^ci = and at tQci = 5. In 
both the CELESTE3D and NPIC simulations, an electron temperature anisotropy develops 
by tQci = 5 both at the center and on the flanks of the current sheet. The anisotropy ranges 
between 2.3 in NPIC and 1.7 in CELESTE3D. The ratio between the ion perpendicular and 
parallel temperature, Tj_|_/Tj||, at tQci = and at tQci = 5, is also shown in Fig. 4. According 
to both PIC codes, ion anisotropic heating is smaller than electron heating, and it is located 
mostly on the flanks of the current sheet. As shown in Fig. 5, the electron distribution 
function is not gyrotropic. In particular, the peak electron perpendicular temperature in 
the z direction (cross sheet direction) is about 25% higher than the temperature in the 
y direction (current aligned direction). Moreover, the heating in the z direction is more 
focused near the center of the current sheet than the heating in the y direction. 

The LHDI also alters the current profile. It thins the current sheet and increases the 
value of the peak current density. The alteration is primarily due to electron acceleration: 
the ion and the electron density modifications are small HHQ, Q- 

In Fig. 6, the 

y-averaged current sheet profile is shown at t = and at tQci = 4 as a function of z. The 
current profiles obtained from both PIC codes agree well and show that the peak current is 
enhanced by about 40% over the initial current profile. 

The LHDI also creates velocity shear. This latter consequence has been the subject of a 
number of papers (see Ref. 4^ and references therein). Its main effect is to promote the 
growth of a KHI, but the three-dimensional simulation shows that the KHI develops after 
the onset of reconnection and thus appears to have little effect on reconnection onset. 



C. Effects of the LHDI on tearing 

Both anisotropic heating and current peaking enhance the linear growth rate of the 
tearing instability and change dramatically its non-linear evolution. The effects of the two 
mechanisms on tearing are considered independently in the following subsections. 
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1. Electron anisotropy 



Figure 7 shows the value of the tearing growth rate 7 as a function of for Tex/Tgn = 
0.9, 1, 1.5,2, and same GEM challenge parameters (in particular Ti/Te± = 5). The growth 
rate of the tearing instability increases appreciably with the ratio Te±/Tf.\\: for Te±/Te\\ = 2 
the maximum growth rate, 7 ~ 3.8^^4, is more than one order of magnitude larger than 
the maximum growth rate for Te±/Te\\ = 1 (7 O.lSQci)- Moreover, the tearing mode 
with maximum growth rate shifts to shorter wavelengths with increasing anisotropy: in the 
case Tf.±/Te\\ = 2 the maximum occurs for fc^L = 2.25, considerably larger that the typical 
kxL = 0.5 seen for isotropic electrons. It should be remarked that electron heating due to 
LHDI does not lead to a simple electron distribution function (see Fig. 5), but assuming 
a bi-Maxwellian distribution function, allows us to estimate the growth rate of the tearing 
instability. 

The tearing eigenmodes for Te±/Tf,\\ = 2 and Te±/T(.\\ = 1 are compared in Fig. 8. 
The most remarkable difference is in the structure of the potential 0, which is larger in the 
presence of anisotropic electrons and with opposite sign. Moreover, for anisotropic electrons, 
the tearing eigenmode involves a narrower region near the center of the current sheet and 
the component is more important. 



These results are consistent wit 
predicted in previous work 



4j,i4a 



1 the enhancement of the tearing instability growth rates 



out our results are obtained by integration 
As an aside, it should be remarked that our linear 



HQ 

over exact orbits, similarly to Ref. 
analysis predicts that ion anisotropy plays only a minor role in thin current sheets (p, ~ L) 
in agreement with earlier work by Burkhart and Chen js^. 

The nonlinear evolution of the tearing mode changes dramatically when the electron tem- 
perature is anisotropic. An initially favorable temperature anisotropy neutralizes the effect 
of the unfavorable anisotropy created by the growth of tearing, which would otherwise cause 
nonlinear stabilization. The evolution of the current sheet is followed in two-dimensional PIC 
simulations in the (x, z) plane, starting with an anisotropic electron temperature. The GEM 
challenge parameters are used with Ti/T(.\\ = 5 on the same domain L^x = 25.6L x 12. 8L, 
but with various electron temperature anisotropies, Te±/Te\\ = 1,1.5,2. In Fig. 9, the 
evolution of the tearing instabilities with mode numbers rrix = 1 — 11 (corresponding to 
kxL = 0.25 — 3) are shown as a function of time. In comparison with the case for initially 
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isotropic electrons (see Fig. 1), high mode numbers grow, forming a number of small tiny 
islands. Subsequently, these small islands merge to form a single large island (the mode 
number = 1 dominates), which encompass the whole domain by tQ^i ~ 20. 

In Fig. 10, the reconnected flux is plotted as a function of time for the three differ- 
ent anisotropies, Te±/T(.\\ = 1,1.5,2, for both NPIC and CELESTE3D simulations. With 
Te\\/Te± = 1, tearing saturates at a low level. For anisotropic distributions, Te\\/Tej_ = 1.5, 2, 
reconnection involves the whole domain and the growth of the tearing instability saturates 
when all the available magnetic flux is reconnected, at a level very similar to the GEM chal- 
lenge. The fast reconnection phase is delayed longer for Tf.±/Tf.\\ = 1.5 than for Tex/Tgn = 2. 

Figure 11 shows a plot of the temperature anisotropy as a function of time for the three 
simulations with different initial anisotropies {Te±/Tf.\\ = 1,1.5,2). The temperature ratio 
{T^yy + Tezz)/'2Texx that corresponds to Te±/Te\\ at t = is plotted. For an initially isotropic 
distribution, Tex/TgH = 1, the ratio decreases as a function of time and stabilizes the tearing 
instability. For initially anisotropic distributions, the ratio Tex/Tgn decreases rapidly to a 
minimum value, Tex/TgH ^ 0.8, and then increases with the onset of fast reconnection. 



2. Current thinning and peaking 

The thinning of the current sheet and consequent increase in maximum current density 
also enhances the linear growth rate of the tearing instability and changes its non-linear 
evolution. This is in agreement with the results of other recent simulation studies 

Figure 12 shows the tearing instability linear growth rate, 7, for different values of the 
current sheet thickness {pi/L = 0.914,1.828,3.656). The growth rate is higher for thinner 
current sheets: 7 = 0.176Qci for Pi/L = 1.828 (GEM challenge thickness) compared with 
7 = 0.632f2ci for Pi/L = 3.656. The maximum growth rate is located at k^L = 0.5 in all 
cases. 

The non-linear evolution of the tearing instability is simulated with both NPIC and 
CELESTE3D. For sufficiently thin current sheets, reconnection is not blocked by non-linear 
saturation of the tearing mode at low levels and encompasses the whole domain. In Fig. 
13, we show the evolution of the current sheet by plotting the reconnected flux, A\E', for 
the Harris sheet with half the current sheet thickness as the GEM challenge {pi/L = 3.656) 
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from the simulations. A domain size x Lz = 12.8di x 6Adi (corresponding to x Lz = 
51. 2L X 25. 6L) is considered. Unlike the standard GEM challenge case without perturbation, 
reconnection involves the whole domain by tQci ~ 20 for NPIC and by tQci ~ 13 for 
CELESTE3D. It should be remarked that the fast reconnection phase is very similar in the 
two codes. The delayed start of the fast reconnection phase in NPIC is probably because 
more particles are used in NPIC than in CELESTE3D, which reduces the noisiness of the 
initial conditions 



IV. CONCLUSION AND DISCUSSION 



Reconnection onset is studied using results from two codes employing very different algo- 
rithms: NPIC is a massively parallel explicit code and CELESTE3D is an implicit-moment 
method PIC code. The results from NPIC and CELESTE3D complement and confirm each 
other. This degree of cross-checking between codes, for which the GEM challenge is de- 
signed, is unusual. Because of their differences in resolution, the agreement between explicit 
and implicit results in detail suggests that the important physical length scales involved in 
the onset of reconnection are comparable to or greater than the electron scales. 

The simulations have pointed out the important role of the LHDI in reconnection on- 
set. Our results confirm that the tearing instability saturates at a low level if no initial 
perturbation is added to the Harris equilibrium. Linear growth of the tearing instability is 
limited by increase in the parallel electron temperature, such that an anisotropy develops 
with Te±/Te\\ < 1 that strongly stabilizes the tearing mode. In three-dimensions, the tearing 
instability does not saturate at small amplitudes. Because the mode spectrum reveals no 
significant oblique mode growth, the current sheet dynamics can be analyzed by studying 
the interrelationship between the current aligned instability that develops in the {y, z) plane 
and the dynamics of the tearing modes in the (x, z) plane. The analysis shows that the 
LHDI strongly modifies the linear and non-linear evolution of the tearing instability. The 
LHDI causes a favorable electron temperature anisotropy, Tex/TeH > 1, and thins and peaks 
the current sheet. Both effects seem effective in enhancing the linear growth rate of the 
tearing instability. With a favorable anisotropy and a thin enough current sheet, the tear- 
ing instability grows large enough to decouple electrons and ions so that reconnection can 
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encompass the whole domain. 

In the current sheet considered there is no background plasma, the guide field is not 
present and no normal component of the magnetic field is introduced. The infiuence of 
these on reconnection onset needs further investigation. 



Both laboratory and satellite observations (32, |33|, |3J] have pointed out the presence of 
the LHDI near reconnection sites, but they have not identified a connection between this 
instability and reconnection through a contributions to anomalous resistivity. 

Satellites have measured electron anisotropy during magnetic substorms. Shinohara et al. 
3^ observe electron anisotropy (Tgx/Ten < 1) at the substorm onset. A more comprehensive 
study of the electron anisotropy during reconnection has been performed by Birn et al. 
IgsI ] . Through one-year average data of satellite measurements, it is shown that an electron 
anistropy (Te^/Ten > 1) precedes substorm onset. At onset, Tej_/Te|| < 1 is observed. 
After onset, Te_|_/Te|| grows again. This behavior recalls closely what is observed in tearing 
simulations (see Fig. 11). Birn a/. ^ also show that an ion temperature anisotropy is less 
relevant than an electron temperature anisotropy, and this also agrees with the simulations 
presented here. 

It is not obvious how to verify current sheet thinning from observations. Our simulation 
starts with a plain Harris sheet from which the LHDI develops and thins the current sheet. 
In reality, on the other hand, the LHDI is always present in current sheets with an amplitude 
that corresponds to the saturation level for that particular current sheet. In particular, MRX 



results j32l. I33| show a thinning of the current sheet, and a consequent enhancement of the 
amplitude of the LHDI instability. It is hard to say if the enhancement of the LHDI is a 
cause or a consequence of the thinning process. 
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Fig. 1: The amplitude of tiie tearing mode are shown as a function of time for the mode 
numbers ttIx = 1 (dashed), rrix = 2 (dash-dotted), rrix = 3 (dotted), and ttIx = 4 (sohd) 
(corresponding to k^L — 0.25, k^L — 0.5, k^L — 0.75, and k^L — 1, respectively), 
during the two-dimensional simulation in the tearing, {x,z), plane. The results are 
from an NPIC simulation. 

Fig. 2: The amplitude of the Fourier modes of the component of the magnetic field 
is shown as a function of the mode numbers m-r and rUy at different times: tfld — 0.2 
(a), tQci = 2.5 (b), tQci = 5 (c), tQci = 10 (d), tQci = 15 (e), and tQci = 20 (f). The 
three-dimensional simulation performed by CELESTE3D is considered. 

Fig. 3: The amplitude of the Fourier modes of the component of the magnetic field 
is shown as a function of the mode numbers m-^ and rUy at different times: tflci — 0.2 
(a), tflci = 2.5 (b), tflci = 5 (c), tfld = 10 (d), tfld = 15 (e), and tfl^ = 20 (f). The 
three-dimensional simulation is performed by CELESTE3D. 

Fig. 4: The enhancement of the electron temperature ratio Te±/Te\\ (a,c), and of the 
ion temperature ratio Ti±/Ti\\ (b,d) due to LHDI is shown from simulations in the {y, z) 
plane. The temperature ratio is shown at time t = (solid line, isotropic distribution) 
and tflci = 5 (dashed line). The results are from NPIC (a,b) and from CELESTE3D 
(c,d). 

Fig. 5: The electron pressure ratios Pezz/Pexx (a) and Peyy/Pexx (b) averaged along 
y are shown at time tfld = 5. The results are from an NPIC simulation in the (y, z) 
plane. 

Fig. 6: The peaking of the current density Jy is shown from the simulation in the 
(y, z) plane. The dotted line represents the current profile at t = averaged along y, 
the profiles at tfld = 4 is shown by the sohd line (NPIC simulation) and dashed hne 
(CELESTE3D simulation). The current is normalized in order that the maximum is 
equal to 1 at i = 0. 

Fig. 7: The growth rate ^/fld of the tearing mode is plotted as a function of k^L for 
Te±/T^\\ = 0.9 (dotted), Tex/Teii = 1 (solid), T^jT^w = 1.5 (dashed), and T^jT^w = 2 
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(dash-dotted). The other plasma parameters are the same as the parameters described 
in Sect. II (in particular, Tj/Tej, = 5). 

Fig. 8: An eigenmode solution for the fastest tearing instabihty is shown. The real 

(sohd) and imaginary (dotted) parts are plotted for (a,e), (b,f), Ay (c,g), and 
Az (d,h) with Te±/Te\\ = 1, k^L = 0.5, and 7 = 0.176flci (a)b,c,d) , and Tg±/Tg\\ = 2, 
k^L = 2.25, and 7 = 3.7820ci(e,f,g,h). 

Fig. 9: The amplitude of the tearing mode is shown as a function of time for the 
mode numbers rrix = 1 and kxL = 0.25 (solid blue), rrix = 2 and kxL = 0.5 (solid 
green), rrix — 3 and k^L — 0.75, (solid red), nix — 4 and k^L — 1, (solid cyan), 
777.3; — 5 and k^L — 1.25, (solid magenta), m-r = 6 and k^L — 1.5, (solid black), 777-r = 7 
and kxL — 1.75 (dashed blue), nix — 8 and kxL — 2 (dashed green), nix — 9 and 
kxL — 2.25, (dashed red), ?77a; — 2.5 and kxL = 10, (dashed cyan), ?77a; — 11 and 
kxL = 2.75, (dashed magenta), and = 12 and k^L = 3, (dashed black), during the 
two-dimensional simulation in the tearing plane with plasma parameters described in 
Sect. II, but Tex/Teii — 2 and Ti/Tg± — 2.5. The results are from an NPIC simulation. 

Fig. 10: The reconnected flux, A^, is shown as a function of time, for simulations with 
the parameters described in Sect. II (solid line), parameters described in Sect. II but 
Te±/Te\\ — 1-5 (dotted), and T^/T^^^ — 5 (dashed), and plasma parameters described in 
Sect. II, but Te±/Te\\ = 2 and Ti/Te\\ = 5. Both NPIC results (a) and CELESTE3D 
results (b) are shown. 

Fig. 11: The electron temperature ratio {Teyy + Tf.zz)/'^Texx (corresponding to Tex/TeH 
at i = 0) is plotted as a function of time, for simulations with the parameters described 
in Sect. II (solid line), plasma parameters described in Sect. II, but Te^_/Te\\ — 1.5 
(dotted), and Ti/Te\\ = 5 (dashed), and parameters described in Sect. II but Te±/Te\\ = 
2 and T^/Tg^w = 5. The results are from NPIC. 

Fig. 12: The growth rate 7/ilci of the tearing mode as a function of kxL is plotted for 
plasma parameters described in Sect. II and different values of Pi/L. 

Fig. 13: The reconnected flux as a function of time is plotted for the simulation in 
the {x, z) plane for a plasma with the plasma parameters described in Sect. II, but 
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Pi/L = 3.656. NPIC (solid) and CELESTESd (dashed) results are plotted. 
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Table I. Simulation parameters 



grid particles AtQo 
(x, z) plane (NPIC) 1280 x 640 16 ■ 10^ 0.03 

{x, z) plane (CELESTE3D) 64 x 64 5 • 10^ 0.45 
(y, z) plane (NPIC) 1600 x 640 2 • 10^ 0.03 

(y, z) plane (CELESTE3D) 128 x 64 1 • 10^ 0.7 
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